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ABSTRACT 

Using results of DNS in the case of two-dimensional homogeneous isotropic 
flows, we first analyze in detail the behavior of the small and large scales of 
Kolmogorov like flows at moderate Reynolds numbers. We derive several 
estimates on the time variations of the small eddies and the nonlinear in- 
teraction terms; those terms play the role of the Reynolds stress tensor in 
the case of LES. Since the time step of a numerical scheme is determined as 
a function of the energy-containing eddies of the flow, the variations of the 
small scales and of the nonlinear interaction terms over one iteration can be- 
come negligible by comparison with the accuracy of the computation. Based 
on this remark, we propose a multilevel scheme which treats differently the 
small and the large eddies. Using mathematical developments, we derive esti- 
mates of all the parameters involved in the algorithm, which then becomes a 
completely self-adaptive procedure. Finally, we perform realistic, simulations 
of (Kolmorov like) flows over several eddy-turnover times. The results are 
analyzed in detail and a parametric study of the nonlinear Galerkin method 
is performed. 


1 Part of the work was done while the second and third authors were visiting the Institute 
for Computer Applications in Science and Engineering (ICASE), NASA Langley Research 
Center, Hampton, VA 23681. 
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Introduction 


A turbulent flow can be characterized by a spatial and temporal chaotic be- 
haviour. Indeed, strong gradients may appear and, in the physical space, 
several structures with various sizes are generated by the external force and 
by the flow itself. The structures are convected by the mean flow. In the case 
of viscous incompressible flows in dimension two, thin and elongated sheets of 
vorticity appear. These structures are characteristic of incompressible flows. 
The phenomenological theory of turbulence, introduced by Kolmogorov in 
dimension three and Kraichnan in dimension two, predicts that the size of 
the smallest scales of a flow decreases while the nondimensional Reynolds 
number grows. So, the number of degrees of freedom required to describe 
the motion can be estimated in terms of the Reynolds number (~ cRe in 
space dimension two). 

Hence, a Direct Numerical Simulation, i.e. the resolution of all physically 
relevant scales, can not be achieved when the Reynolds number becomes too 
high and then, when the turbulence is fully developed. A model is then used 
in many simulations ; i.e. the small scales are not exactly integrated but 
their interaction with the large ones are taken into account in a simplified 
way. Basically, the aim of these models is to recover the large eddies of the 
flow (or their statistic) without explicitly computing all the scales of motion. 

In relation with recent developments in the theory of dynamical systems 
and its application to turbulence phenomena, new objects have been in- 
troduced (exact or approximate inertial manifolds, see Foias, Manley and 
Temam [1]) and new numerical methods have been proposed (the nonlinear 
Galerkin methods, see Marion and Temam [2] and [3], the incremental un- 
known method, see Temam [4]). These objects and methods are based on a 
decomposition of the unknown, such as the velocity field, into its small scale 
component and its large scale one, and essential in the nonlinear Galerkin 
method is a systematically differentiated treatment of the small and large 
scales. 

These approximate inertial manifolds are subsets of the phase space and 
consist in an approximate form of the small scale equations. They provide 
an adiabatic law modeling the interaction between the low and the high fre- 
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quency components of the flow ; the small eddies are in fact expressed as 
a nonlinear function of the large ones, they are slaved by the large ones. 
Moreover, these manifolds enjoy the property that they attract all the orbits 
exponentially fast in time and that they contain the attractor in a thin neigh- 
borhood. In that sense, they provide a good way to approximate the solutions 
of the Navier-Stokes equations. The nonlinear Galerkin method, proposed 
by Marion and Temam [2], consists in looking for a solution lying on these 
specific subsets of the phase space. Several implementations of this scheme 
have previously been done by Jauberteau, Rosier and Temam ([5] and [6]), 
and Dubois, Jauberteau and Temam ([7]). Results presented in these papers 
were mainly feasibility tests involving exact solutions (analytical), where the 
authors tried to recover known velocity and pressure fields. These tests pro- 
vided good indications on computational feasibility and performances, but 
had no physical relevance. In the present article we intend to develop further 
the study of the numerical implementation of the nonlinear Galerkin method 
for flow problems. We have several objectives which we describe hereafter 
in some details. Firstly, we would like to compute physically more relevant 
flows, namely here Kolmogorov type flows. Secondly, the effective implemen- 
tation of the nonlinear Galerkin method involves several parameters (such 
as the cut-off wavelength between low and high frequencies, and the time 
during which the high frequencies are allowed to be frozen) ; and another 
of our aims is to conduct a parametric, study of the method, and to develop 
simple self-adaptive methods for the determination of these parameters. 

Here, we consider two-dimensional periodic flows governed by the incom- 
pressible Navier-Stokes equations. Of course, such flows correspond more to 
an idealized situation rather than a physical one. Nevertheless, this model 
contains several difficulties which occur when simulating turbulent flows, and 
then it provides a good test to implement a numerical algorithm. More phys- 
ically relevant situations, such as three-dimensional flows, will be presented 
elsewhere. In a parallel effort, we are also treating the case of a flow driven 
by a constant pressure gradient between two infinite parallel plates, namely 
the channel flow ( see for instance [8]). 

In the first section of the paper, we introduce the equations and several nota- 
tions. Then, we define the decomposition of the velocity field into large scale 
components and small scale structures. This decomposition depends on a 
cut-off value, corresponding to a wave-number of the Fourier decomposition 
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of the unknowns, and provides two grids in the physical space. In previous 
theoretical works (see for instance [1] and [9]), the authors have rigorously 
proved that, for a sufficiently large wavenumber, the large scale components 
contain most of the energy of the flow. The small scale structures repre- 
sent then small quantities, which have to be taken into account in order to 
correctly compute the large scales of motion. After recalling these results 
and comparing them with numerical simulations, we show, based on both 
theoretical and numerical investigations, that the time variation of the small 
eddies over one time step can be smaller than the energy carried by these 
scales themselves, when the cut-off value is chosen sufficiently large. A sim- 
ilar result is also proved for the nonlinear interaction terms expressing the 
action of the small eddies over the large ones. Then, introducing the expected 
accuracy of the computation £ as a parameter, we deduce that there exists 
a level such that the one-step time variation of the small scale components 
is smaller than £. We note here that this result is not in contradiction with 
the fact that the small eddies evolve faster than the large ones. In fact, their 
time scales are smaller but their order of magnitude also decreases for large 
wavenumbers. 

Hence, the time variations of part of the spectrum can be locally neglected. 
We have also proved that the interaction terms enjoy the same property. 
These two fundamental results are one of the keys of the numerical algo- 
rithm proposed in the following sections. 

Based on these results, we have implemented a spatial and temporal mul- 
tilevel adaptative method ; i.e. the level of refinement, which define the 
separation of the flow into low and high frequencies, evolves in time under- 
going successive multilevel V-cycles ; moreover, it follows the dynamic of the 
small and the large scales of the motion. Moreover, some part of the spec- 
trum are frozen, as well as the interaction terms, during short time periods. 
After a complete description of this multilevel procedure, we derive several 
estimates of characteristic time scales for the small structures and for the 
terms involving their interactions with the large scale components. These 
estimates provide an efficient way to evaluate the length of the time interval 
during which the small scales, corresponding to several levels of refinement, 
as well as the corresponding interaction terms are allowed to be frozen with- 
out introducing an error larger than the accuracy e. Hence, we obtain a 
completely self-adaptative procedure which depend on only one parameter, 
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namely e. At the end of each frozen period, the velocity field is projected 
on an Approximate Inertial Manifold (AIM), which provides a simple and 
efficient way to directly compute the high frequencies. We analyze here the 
perturbation introduced by using first order Approximate Inertial Manifolds 
and we deduce in which range of the spectrum the manifold can be used to 
compute the small scales. 

The last section is devoted to the description and the analysis of several 
numerical simulations performed with the multiscale (Nonlinear Galerkin) 
method. The integral scale Reynolds number Rei, was successively set to 
784 and to 6,328 for spatial resolutions of respectively (256) 2 and (5 1 2 ) 2 . 
In these cases, a full dissipation range is resolved as the larger wavenumber 
is at least 6 times larger than the dissipation wavenumber. The Reynolds 
numbers are not sufficiently large for the solutions to have a k~ 3 energy 
spectrum power range ; indeed, the computed velocities have an intermedi- 
ate energy spectrum between k~ 4 and k~ 3 . Nevertheless, the different scale 
components are strongly time dependent, and they provide interesting tests 
for the multilevel adaptative procedure. More physically relevant situations, 
eg simulations with larger Reynolds number, will be treated and presented 
elsewhere. For each simulation, the code has run during more than 50,000 
time iterations, which represents 10s of unit hours on a Cray-2. 

In this last section, we first compare the results obtained with a previous 
version of the algorithm with those obtained with the multiscale procedure 
described in subsection 2.4. For this comparison a Direct Numerical Sim- 
ulation performed with a standard pseudo-spectral Galerkin approximation 
is used as a reference. In the previous versions of the algorithm ( [7], [10]), 
there was no control of the variations of the small eddies. Hence, the cut-off 
wavelengths as well as the length of the period during which the variations 
of the small scales are frozen were not properly evaluated. This led to an 
accumulation of perturbations and a loss of the accuracy ; this is no more 
the case. 

Finally, in order to study the effect of the cut-off wavelength on the numeri- 
cal results, we perform several simulations with decreasing cut-off values. In 
fact, we show that the separation wavenumber can be taken of the order of 
the dissipation wavenumber k v ; as k v is of the order of 20, in our computa- 
tion, the cut-off value can not reasonably take smaller value. Nevertheless, 
the standard Galerkin pseudo-spectral method, with a larger wavelength of 
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the order of the dissipation one, does not provide a sufficient resolution : the 
large scales of motion are not recovered in that case. Hence, the multiscale 
method provides an efficient way to simulate the enstrophy dissipation and 
the energy dissipation ranges. In order to describe the effect of the multiscale 
procedure inside the power range, simulations with larger Reynolds number 
have to be performed. Such study will be presented elsewhere. 
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1 Motivations. 

1.1 Preliminary. 

We consider here two-dimensional viscous incompressible flows driven by an 
external volume force f, and governed by the Navier-Stokes equations : 

Qjf - I'Au + (u • V)u + Vp = f, 

< Vu = 0 , ( 1 ) 

. u(x, t = 0) = u 0 (x), 

where u(x, t) = (ui(x, t), u 2 (x, t)) is the velocity field, p(x, t) the pressure 
(x = (xi,x 2 )), v is the kinetic viscosity and the density is set to unity. 
Equation (1) is supplemented with boundary conditions, namely u and p are 
periodic of period L j (resp. L 2 ) in the direction X] (resp. x 2 ). We denote 
by ft = (0, L\) x (0, L 2 ) the period and we assume that the average of the 
external force over the period vanishes, i.e. : 

— f f(x,t) dx = 0 , Vf. 

I SZ | Jn 

Taking the average over the whole domain ft of the momentum equation 
in (1) and using the periodicity of u and p, we obtain : 



Assuming now that the average of the initial condition is also zero, we con- 
clude that the average of the velocity field vanishes : 

7777 / u(x, t) dx = 0, Vf. 

\ ll \ Jn 

We now denote by uj = V x u the vorticity and as usual we rewrite the 
conservation of momentum equation in (1) as : 

du 

— -vAu + (ux\i ) + VP = f, (2) 
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where P = p + ^ | u | 2 • This choice is motivated by two facts. First, it 
is important to note that a pseudo-spectral evaluation of (a? x u) requires 
three Fast Fourier Transforms (F FT) fewer than the usual form (u • V)u. 
Also, this form of the equations semi-conserves the kinetic energy for inviscid 
flows, when a collocation or pseudo-spectral discretization is used (see [11]) 
and ensures numerical stability, while the original form does not. 


Since the unknowns u and p are space periodic functions, they can be 
expanded in Fourier series : 


u(x,/) = J2 u(k,0 lt, k(x), 

keZ 2 

p(x,<) = p(M) ™k(x), 

keZ 2 


(3) 


where w k (x) = c tarkt-x , k L = ( fc, / L\ , k 2 /L 2 ) and k L x = {k^x^/L^ + k 2 x 2 /L 2 ) • 
The Fourier coefficients u(k ,t) = (ui{k u t),u 2 {k 2 ,t)) (resp. p(k,t)) are com- 
plex numbers and we recall that : 


= u,(k,i), i = 1 or 2, 
P(~k,t) = p(k ,<), 


(4) 


where c denotes the complex conjugate of c. This property is very important 
in practice ; it allows to store twice fewer coefficients in the second direction 
of the spectral space, when we look for a finite approximation of u(x, t ) and 
p(x, t), i.e. by storing u(k,<) for k such that k 2 > 0, (4) allows us to recover 
the Fourier coefficients u(k,t) for k 2 < 0. 


Let us now introduce Pdiv the projection operator onto the divergence 
free space ; Pdiv can be easily expressed as : 

PdivV>(x) = £ (<£k- TTTl( k -^)] ^k( x )- 

kgZ 3 t I k 1 ' 

Assuming that u and p lie in the proper Hilbert spaces (see for instance [12], [13]) 
and applying Pdj v to the Navier-Stokes equations, we obtain 

^ - i^Au + P(u,u) = g, (5) 

at. 
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where g = Paj v f and f?(u,u) is a bilinear form defined by : 

B{ u,u) = P dlv (u> x u) 

= E (( WX U )k - 7~j~f2 (k • lU k (x) 

keZ 2 V 1 k 1 ' 

As we shall see later, the numerical procedures are directly applied to this 
last form of the Navier-Stokes equations. 


1.2 Separation of scales. 

1.2.1 Small and large scales equations. 

Let us choose an integer N. We introduce fV the orthogonal projector onto 
the space spanned by the first N 2 Fourier modes. A pseudo-spectral Galerkin 
approximation Ujv of the velocity field u is given by 

u/v(x, t) = J2 u(k ,t)w k (x), 

ke/ N 

where I N — [\ — N/2, N/2] x [0, A^/2] , and satisfies the following equation 
dujv 

— vAu n + P n B(u n ,u n ) = P N g. (7) 

Here we only consider external force with no high wavenumber coefficients, 
namely : 

g(k) = 0, for k € I N \I mo , for some m 0 . (8) 

Instead of (8), we can also assume that the spectrum of the external force 
has an exponential decay for wavenumbers larger than a given one. 

We now introduce the decomposition 

Uat = y Nj + zyv, , with m 0 < TV, < N, (9) 

where 

y/v, = Pn^n = 53 u(M)u> k (10) 

k 6/^ 
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and 

Zjv, = Qn^n = 53 “(MW- (11) 

k€/^=/Ar\/jv, 

We note that, in the decomposition (9), yjv, represents the large-scale struc- 
tures of the flow, and corresponds to the small-scale components, or in 
other words, the fine structures of the flow. Since, the operators P/v, and 
Qjv, commute with the differentiation operators, we have : 

Pjv,( Au/v) = A(P Nl u N ) = Ayjv,, 


and 

= A(<5jv, u /v) — Azjv,. 

By projecting (7) with respect to Pat, and Q/y,, we obtain the following 
coupled system of ordinary differential equations (ODE) : 

- vAy Nl + PN!B(y Nl +z Nl ,y Nl +z Nl ) = ^V,g, 

< (12) 

+ Q^ i B(y Nl + z WiJM + zyv.) = 0. 

The initial conditions associated with (12) are P/^Uo and Qn, u o- l* 1 (12), as 
well as in the following, we will only use the projection operator Pjv, with 
N\ larger than mo- 

We now introduce two norms : 

I I2 = (J a I <p(x,t) | 2 dx) , 
and 1/2 

IMI = (jf | V V (x,() | ! <Jx) \ 

for any given field <p(x, t) = (^i(x, t ), ^(x, t)). We note that these norms are 
related to well-known physical quantities. Indeed, with the above definitions, 
the kinetic energy e(u) of the velocity field u is related to | — (2 by 

e(u) = l - | u \l, 
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and in the case considered here of periodic boundary conditions, the enstro- 
phy, en.s(u) is : 

U It = \ I - \i . 

We now recall theoretical results, constituting the basis of this framework ; 
they were established by Foias, Manley and Temam in [1] and [9]. In [1], the 
authors proved that for sufficiently large values of N and N\, namely N and 

N\ ~ G’ 2 , where G designates the Grashoff number G = ^ f J 2 (A] ~ l/L 2 

V A\ 

where L is a characteristic length of the domain), and after a transient time, 
depending only on the data, the small scale components Z/Vj (t) remain small 
in both norms introduced above, while the large scales y/Vj(0 are of the order 
of u (t). We do not know how close to the inertial range the corresponding 
cut-off eigenvalue \n x — cN 2 can be, and it is one of our objectives in this 
work to explore this point. In [9], the authors derived other estimates of the 

k d 11^,011 
|ywi(Ola Ilyjv.WII 


ens(u) = -| 


in terms of the Grashoff number at appropriate powers of l/( N\ + 1), when the 
cut-off value N\ is still of the order of G 2 . Moreover, by using the analyticity 
in time of the solutions of the Navier-Stokes equations and invoking the 

Cauchy formula, one can show that | — y | 2 , as well 11 


dt 


as 


w i 

Ht 


are 


respectively of the order of | z^ I 2 and || z ^ ||, for N\ sufficiently large, 
corresponding to a wavenumber in the dissipation range. Even if these results 
do not provide fine estimates on the norms of the time derivative of the small- 
scale components, they show that the time derivatives of the velocity field u 
have also a decaying spectrum, which is not the case of course in the inertial 
range. 

This behavior of the small-scale z^ is very important and is one of the 
keys of the numerical implementation of the multilevel method, as we will 
see in the following sections. Although, these results were proved theoreti- 
cally only when is of the order of G 2 , and is far in the dissipation 
range, numerical results show that Zyvj is small compared to for much 
smaller values of , and the same is true for their time derivative. From 
the phenomenological theory of two-dimensional turbulence point of view, 
due in part to Kraichman (cf [14] and [15]), the energy spectrum of a de- 
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veloped turbulent flow is expected to display a power-law inertial range and 
an exponential decay for wavenumbers larger than the Kraichnan's dissipa- 
tion wavenumber (see Figure 1). So, the theoretical results mentioned above 
are in agreement with this assumption in the sense that for sufficiently large 
wavenumbers, the energy spectrum has a decaying behaviour. The exponen- 
tial range is generally referred as the dissipation range. 



Figure 1: Energy Spectrum E(k). 

Due to the bilinearity of B, we can split B( Ujv,U/v) into 
B( U/v,U/v) = 5(yNi + z Ni i y^i + Z N! ) 


— 5(y/V,,yjV,) T fiint(yA/i ? Z N] )? 


(13) 


where 


fiintiywpZiv,) = B{ y/VpZ/Vi) + B{z Nl ,y Nl ) + B{z Nl ,z Nl ). 

It represents the interaction between the small and the large structures and 
the interaction of the small structures among them. We now recall that a 
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classical (Galerkin) approximation of u based on the first iV 2 Fourier coef- 
ficients consists in setting z N] to zero in (12), so that P Nl B mt (y Nl , z Nl ) is 
neglected. If zyv, represents a small part of the kinetic energy, namely 

e( z /Vj ) < e(y Nl ), 

then, we can expect that 

I F , N 1 fii„t(yN 1 , z Ni) I 2 < I FW, B(y Nl , y Nj ) | 2 . 

Therefore, the contribution of the interaction terms is much smaller than the 
component of the nonlinear terms, involving only the large scales yjv,. This 
argument can be justified when comparing Figures 4 and 5. More precisely, 
we have the following 

I I2 < | Pn, B(y Nl ,y Nl ) | 2 


v | Ayw, | 2 


Nx 

dt 


However, the interactions terms can not be neglected in the first equation ( 12) 
unless N x is taken very large. Indeed, their effects on long time computations 
can modify the behavior of the large scale components. Since the evaluation 
of these terms at each time step is very expansive, we try to take their effects 
into account in a simplified manner. We will see that the time variations 
of Zjv, and FV, F?j nt (y/v, , zjv, ) are locally negligible, and this will allow us to 
freeze them during some parts of the iterations. 


1.2.2 Time variations of z^, and FV, ^^(yiVnZivJ- 

We now aim to derive an estimate of the variations of z^r, over one time 
iteration. This quantity can be represented by 

A W] = At | zjv, | 2 , 

where the dot represents the differentiation with respect to t. The time step 
At is given by the CFL stability condition : 

At N | u;v |l°° < <*, with a < 1. (14) 
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We recall that N is the total number of modes in each direction. We assume 
that the smallest scale (n = l/fcyv, k N = -W/2, is smaller than the Kraichnan 
dissipation scale t v — l/k v . Let us first consider the case where = 1/^ is 
lying between i v and i N , and i N , f Nl , U are of the same order of magnitude : 


f-N < ^Ni < 

Let us assume, as shown in Figure 6, that the time derivative z^v, is of the 
order of the dissipation term (remember we are in the dissipation range), 


| ZN, b ~ " I Az^ b • 

Due to the exponential decay of the velocity spectrum in the dissipation 
range, we can write 

v | Aza?! | 2 < ci vk 2 Nx | zjvj I 2 , (15) 


where Ci is a nondimensional constant of the order of unity. It follows by (14) 
that : 

Aa/j = Af | za/j I 2 5: c 2 v/\t k Ni I zat, I 2 

(16) 

5 c * nj ign ^^' 1 * n ‘ u ■ 

Then, using the estimate of k n obtained in [16] namely 


k n ~ k 0 G 1/3 = 


1 . 1 g ii /3 


we obtain : 


Aa/i < c 2 


Q^O 1 g l2 /3 
v^A 1 / 3 | ua/ U°° 



As < A: at, we obtain : 




^ £3 1 


1 1/3 
12 


Ujv | L°° 


g 


©)• 


1/3 I ZNi I 2 


Since k Nl ~ then for sufficiently small values of the viscosity 1 /, the vari- 
ations of over one time iteration are much smaller than | Zjv, I 2 • 
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We now try to derive a similar estimate when 


A; < 

i.e. when f. Nl is in the inertial range. In that case, we assume as shown 
on Figure 7, that the time derivative z Nl is of the order of the interaction 
nonlinear terms 

| Z A r , I2 ~ | Qn 1 7 ?int(y/Vi 1 Zyy, ) I2 • 

We recall that 

I Qn, #int(y/Vi , Z/V, ) I 2 


= I 5(y/v, , Z /V, ) + Qn x B( Zyv, , yyv, ) + Qn, B{z Ni , Z N] ) [2 

(17) 

^ I Qn, B{yN 1 , Zyv ( ) I2 + | y/Vj) I2 

+ I Q/V, ^( z Ah> z Ni) I2 • 

As Qjvj > s a projection operator, we have 

I B{y N , , z/Vi ) I2 < | B( y^,,z/v, ) I2 • 

The bilinear form 5 can be estimated as follows (see for instance [12] and [13]) 

I B{ y/VnZ/v,) | 2 < C4 | y/v, |l°°|I z Nx II • 

As | yty |l°° — | l^oo, for sufficiently large values of AT 1? we obtain 

I B(y Nl ,z Nl ) | 2 < c 4 | | £/°° 1 1 Ztv, II ■ 

The decay of the velocity Fourier components implies that 


I z yv, || < c 5 k Ni | z Nl | 2 , where c 5 ~ 1, 


so that 


We also have 


Qni B{yN i i Z/Vj ) I 2 < c 4 c 5 k Ni | l^oo | Z Nl | 2 . 


QN t B{z Nl ,y Nl ) | 2 < | B(z Nl ,y Nl ) | 2 , 
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which can be estimated by (see [12], [13]) 

| jB(zyv, , y/V, ) 1 2 5; c &\ Z AT, I I Az/v, \2 IlyN.II- 

In fact, the norm || y/v, || and the infinity norm | y yv, |l°° are of the same 
order. We then obtain 

| B(zyv,,yN,) |2 < c 7 | Z N] fi /J | Azyv, ll^l y/V, |l°° - 

From the inequality (15), we deduce that 

| B(zyv, , yyv, ) I2 < C 8 fcyv, | y yv, |l°°| Zyv, I2 • 

A similar estimate can be derived for the third term B(z ^ l , z/v, ). We 
finally obtain 

| Qj^fiintty/VnZN,) I2 < C 9 fcyv, I U N IL-Izyv, I2 • 


Then, 


As 


we find 


Ayv, < c 9 kyv,Af | uyv j^oo | Zyv, I2 


Af < 


a 


a 


UjV U<°° 


AT 


UJV |l 00 \/2 Aryv 


Ayv, < c 9 


a / 1 
\/2 1 &yv 




zyv, I21 


which also implies that Ayv, <| zyv, I2 in the inertial range as well. Fig- 
ure 8 supports the theoretical estimates derived above. We want to note, 
at this point, that these results will remain valid if energy is injected in the 
small scales by a given external force, having a decaying spectrum ; i.e. if 
g(k) ~| k |“ e-M, for | k | larger than a given wavenumber. 


Let us now introduce e the accuracy of the computations ; e represents 
here an energy error. We assume that £ is a given parameter, which can 
be chosen under several considerations ; examples will be given in Section 3 
devoted to the description of the numerical results. The accuracy e can be 
associated with a small scale, i.e. 

there exists a level N\(e) (< N) such that 
| Zyv,( e ) £• 
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According to the estimate previously derived, Ajv,( £ ) is much smaller than e. 
From the previous estimates we can deduce that there exists several levels of 
discretization N\ lower that A^(e), for which A n, is also smaller than e. Let 
us denote by N[(e) the lowest level of refinement satisfying Ajv'( £ ) < £. The 
above result means that the small scales can be integrated with a larger time 
step, even if their time scales are smaller than the ones of the large eddies. 
By definition and due to the behavior of | z/v, | 2 , with respect to , the 
quantity Ajv, increases for decreasing values of N^. Hence, the appropriate 
time step for ^Vj( e ) is larger than the appropriate one for Moreover, 

on Figure 4, we note that | Zyv, | 2 presents strong variations during the time 
evolution. Consequently, the levels N^e) and N[(e) are not constant in time. 
To take into account this specific dynamics of the small scales, we propose, 
as it was done in [10], [7] and [17], to use multigrid technics. The integration 
of the intermediate scales F/v t1 ^Ni(e) < f^/Vi < Cvj ( £ ) , consists in performing 
a succession of multigrid V-cycles during which some parts of the spectrum 
are frozen during the time evolution. Hence, we use here a property of the 
first order term of the evolution equation of the small scales. 

On Figures 1 1, we can see that the variations of P Nl B int (y Ni , z Nl ) and of 
z Nl are correlated. Therefore, we deduce that the variations of P Nl £?i nt (yjv, , Zjv, ) 
over one time iteration is smaller than F/v, Ffin^yw^Zw,) ; in fact, one can 
show that : 

&'ni = Af | P Nl B int (y N] ,z Nl ) | 2 < c 10 I FW, ^(y/v, , z/v, ) | 2 • 

The nonlinear interaction term P Nl B int ( y Nl ,z Nl ) represents a small part of 
the large scales time derivative y ' Ni . Hence, its variation over one time itera- 
tion can be neglected without a lost of the accuracy on the large scale approx- 
imation. We use here a property of the time evolution of P Nl 5i nt (y^, , z^, ) 
which is a second order term of the large scale evolution. As the time scale of 
z^Vj is much smaller than the one of z^, , the time derivative P/v, Pint(y/Vi i Zyv, ) 
behaves as Zyv, and then is very oscilating with the time. From the previous 
inequality, it appears that the quantity A' Ni can be controlled by the term 
FW, B int (y Nl ,z N] ) ; hence, we estimate the level N[(s) by imposing that the 
ratio 

I PNl (t) Pint ( y JV' (e) ) Zyy> ) 1 2 

I Pn[(c)B( y N[(c),y N[( e)) b 
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remains smaller than a given constant. As the computation of the nonlinear 
terms requires several calls to the FFT routines, the evaluation of the above 
ratio is expensive and then prohibitive. So, we have looked for another form 
of it. In the numerical results, we have noted that 


P N\ N \ } ^N\ ) 

2 ||*Wj | 

- 1 

ens(zyv, ) \ 

I Pn x B(y N} ,y Nl ) | ; 

2 IIyn.II ' 

lens(y/v,)/ 


as we can see on Figures 13. Hence, we use the ratio of the enstrophy of the 
small scales over the enstrophy of the large ones to evaluate the level N[(e). 
As we have seen before, the notion of variation of the scales is intrinsic to the 
definition of N[(e). More sophisticated criteria will be derived in Section 2.2. 
Finally, we want to note that controlling the size of P iv, #int(yw, , ZNi ) a °d its 
time variation ensures that the interaction of the small scales over the large 
ones is negligible, in the sense that this interaction can be locally neglected 
from a numerical point of view. 


2 Description of the multiscale method. 

As it was shown in Section 1.2.2, the small scales of the flow as well as the 
terms involving their nonlinear interactions can be fixed in time during a few 
time steps. Nevertheless, their order of magnitude may change drastically 
over a period of time ; so, the cut-off value A^i defining the separation between 
the small and the large eddies can not be let fixed in time. Hence, we 
propose a multilevel adaptative procedure evaluating the appropriate level 
of refinement as time evolves. Using partly theoretical arguments, we show in 
this section that we can derive estimates for the variations of the small scales 
and of the transfer terms to the largest scales. We then deduce estimates 
for the length of the frozen periods. Finally, we introduce these estimates in 
the algorithm and we derive a dynamic procedure allowing an a priori and 
an a posteriori control of the length of the quasi-static time intervals. In the 
first subsection, we describe the multilevel adaptive procedure ; secondly, we 
derive time scales estimates of the fine structures and of the transfer terms. 
In the third part, we give an explicit approximation law used to compute 
the small scale components of the flow and we derive error estimates for this 
Approximate Inertial Manifold. Finally, we describe the whole algorithm 
including the modifications previously discussed. 
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2.1 The multilevel adaptive procedure 

As in the preceeding Section 1.2.1, we choose an integer N larger than m 0 
(cf (8)), which represents the total number of modes retained in the trunca- 
tion and we denote by 

u * = E u(M) w k , 

ke/;v 

an approximation of u. We associate with N the largest wave number of 
U;v, namely kpj = N/2. Hence, the smallest scale in the computation is 
= 1/&N- 

We are now given a sequence of levels Ni satisfying 

N-i < N 2 < ... < Ni < /V 1+1 < ... < N. (19) 

As we want to perform pseudo-spectral approximations of equations (12) on 
these different levels of refinement, the elements N{ of these sequences have 
to satisfy the restrictions imposed by the Fast Fourier Transforms (FFT) ; 
namely, the N^s must be of the form 2 P 3 9 5 r , where p > 2 and q, r > 0. Such 
an algorithm enables us to define a suitable sequence of levels Ni. Examples 
of sequences will be given in Section 3, where the numerical results will be 
described. We note that an FFT allowing only decompositions in powers of 
2 is not efficient for this purpose. 


Let us now assume that the approximation Uyv(x, t) is known at a time 
tj. As it was suggested in the previous section, we define two levels of dis- 
cretization A r q(fj) and N i2 (tj) by the following procedure : 


i\ is defined by the condition that 
for every i > i \ , 


ens(z ; y t (f ,)) 
ens(y Ni (t])) 


< 0i. 


?2 is defined by the condition that 
for every i > 


e (y n,(^)) 


< 02 • 


( 20 ) 


( 21 ) 
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where 6 \ and 0 2 are two given constants ; 62 is chosen so that e(z/v, 2 ) is of the 
order of the accuracy e 2 . In the previous versions of the algorithm, the pa- 
rameter 9\ was arbitrary fixed ; we will derive an estimate of 0 j in section 2.4. 
. In order to be sure that TV,-, < Ni 2 , we impose the additional condition that 
(i -2 ~ * 1 ) has to be larger than a given constant. This is motivated by the fact 
that, as we described in the previous section, the intermediate region of the 
spectrum between A and N t 2 (tj) is a transition zone between the large 
scales and the static (frozen) small scales. 


We now introduce the quantities 


and 


e x (t 3 ) 


/ ens(zA^jt^y /2 

\ens(y Nii {tj))) 


Htj) 


V e (y N, 2 (tj))J 


(22) 


(23) 


Let r c (tj ) be the length of the time interval during which the scales smaller 
than , i.e. Z/v, 2 , can be frozen without loosing the order of approximation 
on the larger scales ; estimates for the characteristic length T c (tj) will be 
derived in Section 2.2. In previous works, see for instance [8], r c (fj) was 
estimated by the characteristic relaxation time of the viscous term, namely : 


Tc(tj) = {vk 2 Ni2 ) l . 

As we will see in Section 3, this estimate is not fine enough and may induce 
strong errors in the approximation of the velocity field. The available levels 
of refinement, lying between Nq (tj) and Ni 2 (tj), are 


JV„ < JV, 1+ i < ... < Ni < ... < N h -i < N i2 , 

which correspond to (*2 — + 1) levels. For the sake of simplicity, we omit 

here the dependence on tj for the levels N li and Ni 2 . 


As in classical multigrid methods, we use the concept of V-cycle to per- 
form the integration of (12) on the interval [tj,tj-\-T c ]. Let us define a V-cycle 
starting at time tj. Such a V-cycle is constituted of two phases described as 
follows : 
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• phase 1 : on the interval [tj,tj + (i 2 — the current level Ni(t) is 

defined by : 

Nil 0 = for j = 0, i -2 — t'j, 

hence Ni(t) decreases from N{ 2 (tj) to TV,-, (<_,). 


• phase 2 : on the interval [tj + (i 2 - ii)At,tj + (2 (i 2 - i'i) + 1)A t], the 
current level N,(t) is defined by : 

Ni(i) = 1), for j = i 2 - i ! + 1 , 2(i 2 - i'i) + 1 , 

hence Ni(t) increases from N tl {t 3 ) to /V, 2 (tj). 


Level N i 



Then, a V-cycle consists in [2(i 2 — »!) + !] time iterations. The quantity T c (tj) 
is adjusted so that the time interval [tj, tj + t c ] can be divided into a fixed 
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number of such cycles. Figure 2 summarizes this process. 

Let t be an intermediate time on the interval [tj,tj + t c ] ; according to 
the previous procedure, the current level N{(t) is given by : 

' N ti . r+ 1 , if l<r<(* 2 -*i), 

Ni(t) = J yv 11+(r _ (t2 _ tl+1)) (24) 

if (^2 — *i + 1) < r — 2(*2 — *l)> 

where r is given by : 

t-tj = (2p (i 2 - fi) + r) At. 

Knowing the size N, of the coarse grid at time t, we decompose u N (t) into : 

Ua l{t) = y Ni(t) + ZN,{t), 

where yN,(t) represents the scales larger than f/v, , and z^(t) the scales 
smaller than t Ni and larger than l N . The computation of both components 
y N .(t) and ZN,{t) are performed as follows : 

• computation of z at, (t) : 

z Ni (i) = z N ,(t-At), 

i.e. za/, (t) is frozen and set to its last value. 

• computation of y : in order to evaluate y Ni(t), we integrate 
equation (12) over the interval [t - At, t] ; then it follows that 

tt(k,f) = e -l 'l k l* A< u(k, i — At) 

+ Jt At e_Hk|2(< ” T)fik ( yNi ( r ) ,yN> ( r )) dT ( 25 ) 

+ f e -‘'l kl2(t - r) B i „ ti k(yA?,(T),z Nl (r)) dr 

Jt-At 

for every k in 7 a/, = [1 — NJ 2, iVj/2] x [0, iV,/2]. 
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The first integral is computed by an explicit Runge-Kutta scheme of order 
3 (see [7]). With this scheme the interval [t - AM] is split into three sub- 
intervals of the form [<„ Mi], where t 0 = t - At and t 3 = t. On each of these 
sub-intervals, the second integral is approximated as follows : 

r +1 e ^ |k|2( ^'- i, J B intik (y Ni (6),z N ,((i)) dS 

J t, 

( 26 ) 

- ,-WWd 

We note here that this approximation is an explicit Euler scheme on the 
interval where the following approximation is performed 

^int,k(yw,(b), ~ fiint,k(y/V,(M, z N,{tj))- 

This integration requires the storage of P Nt B int (y Ni (tj), z N , (tj)), at the be- 
ginning of the cycle, for each coarse grid TV, between TV,, and N, 2 , i.e. for 
(i% — *i + 1) levels. 

As the time T c (tj) is adjusted to be a multiple of a complete V— cycle, 
the current level at time tj + T c(tj) is equal to N t , 2 , the highest coarse level. 
Referring to the integration process described above, the large scales y^v 
are known at the end of the cycle, i.e. at time tj + r c (<j)- The smallest scales 
z n, 2 are then updated by projecting the approximate solution u^(tj + r c (tj)) 
on an approximate form of the small scales equation (12), for instance 

ZN i3 (tj + T c (tj)) = <f>{yN. 2 {t J + T c {t ] )),Z Nti {t J \T c {t J )). (27) 

(27) is the equation of an Approximate Inertial Manifold. Such manifolds 
were first derived in [1] ; (27) provides an interaction law between the small 
and the large scale components of the flow, and expresses z Ni as a function 
of y n, 2 • Further information on such law can be found in [18], [19], and [7], 
Other kinds of approximate inertial manifolds are derived in [20] and [21] ; 
the implementation of an Approximate Inertial Manifolds of first order will 
be discussed in Section 2.3. In Section 2.3, we will discuss on the efficiency 
of these nonlinear forms and derive some estimates which explain in which 
range of the spectrum they can be implemented. From a strictly compu- 
tational point of view, an approximate inertial manifold is efficient in the 
sense that this equation allows us to estimate the small scales as an explicit 
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function of the large ones. Moreover, a first order law only requires one eval- 
uation of the nonlinear terms, on the fine grid (see Section 2.3). 

Finally, at t i+ 1 = t, + r e (tj), i.e. at the end of a whole cycle, we have an 
approximation of ujv(tj+i) by : 

UAf(fj+i) = yjv, 2 (<j+i)+ z iv, 2 ( i j+i)- 

We start the full procedure again by computing two new levels N n (t J+ \) and 
Ni 2 (tj+\). Then, we perform new V-cycles on the time interval [tj+i,tj+i + 
r c (tj+ 1 )], and so on. 

Basically, we can summarize this process by saying that the full time in- 
terval [0, T] of the whole computation is split into several small time intervals 
[tj,tj + T c (tj)\ and, on each of these intervals the velocity spectrum is split 
into three fundamental regions : 

• the dynamical range, corresponding to wavenumbers smaller than , • 

It represents the large scale structures of the flow - i.e. the scales con- 
taining most of the energy and the enstrophy of the flow. The modes 
lying in this part of the spectrum are integrated at each time step of 
the time interval + Tc(ij)]. 

• a transition range , corresponding to wavenumbers between fc/v,, and 
k N . . An up and down oscillation process is used, i.e. the current level 
of discretization Ni(i) undergoes all the intermediate levels between 
N tl {t) and /V, 2 (t), while the time evolves. 

• a quasi-static range, i.e. for wavenumbers larger than k n, 2 - It repre- 
sents the smallest scales, which are numerically negligible ; i.e. their 
energy and their variations are smaller than the expected accuracy. 

Figure 3 represents these three different regions of the spectrum. 

Now we want to make some remarks on the effect of the integration 
procedure previously described on the smale scale components of the flow. 
One of the crucial points concerns the technic used to update the small 
structures z Nl (t) lying in the transition range. 
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Dynamic 



Figure 3: The multilevel procedure. 


Let us first recall the equation satisfied by the k ^ coefficient of the Fourier 
decomposition of Ujv : 



+ v | k | 2 u(k) 


+B\t{yNi,yN l ) 


+ 5i„t,k(yA? 1 ,ZA? 1 ) = o. 

We assume here that k is such that | k |> m 0 and hence (see (8)) : 

g(k) = 0. 

We can then write (28) in the following form 

j t (e-W^ifk)) = -e Hkp ‘ [%,„>,) - Bta«.k(yn,,z n ,)] . 
Integration of (29) over a time interval + r] leads to : 
u(k ,t + r)= e- 1 / ' k l 2 T u(k,<) 

— J t ^ t+r ^ [fik(y»»i 5 y ni ) + Sint,k(ynj , z„, )] do 


(28) 


(29) 


(30) 
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We note that if k is large enough, u(k) is a component of the small scales. 
For instance, we assume that : 

I k |> N h {t). 

Then, k lies either in the transition or quasi-static, range. According to the 
implementation of the multilevel method in (30), u (k,<) is replaced by an 
approximation close to the actual value, as the perturbations are smaller 
than the accuracy e. If k lies far enough in the dissipation range, then the 
dissipative terms are quasi— dominent, namely : 

I Zjv, | 2 I t'Azyv, h - | QNiByniiyNnZNt) U, 

as we can see on Figure 12. In that case, the errors introduced by the quasi- 
static integration are quickly damped by the effects of the operator A . 
few corrected iterations are needed. If k is not as large, i.e. if it is closer to 
the inertial range, then the coupled nonlinear terms become the most impor- 
tant terms of equation (30). In that case, the error is convected by nonlinear 
effects in the largest wavenumbers, and are finally damped by viscous effects. 

During the V-cycles, the modes closer to jV„(f) are more often integrated 
than the ones close to N, 2 (t). The structure of the V-cycles is thus well- 
adapted to the integration of the small scales. As we will see in Section 3, 
devoted to the description of the numerical results, there is no accumulation 
of errors in the intermediate scales and there is no energy pile-up in the high 
wavenumbers of the Fourier decomposition. The enstrophy cascades are well 
described by the V-cycle technic. 


2.2 Time scale estimates for the small eddies and for 
the nonlinear interaction terms. 

Tim e scale estimates for z^. 

From now on and for the sake of simplicity, we omit the subscripts Ni in the 
notations. We then denote by z (t) the small scale components of the flow. 
As the Navier-Stokes equations are analytic in time (see [1]), we can write a 
Taylor expansion of z(t) : 

z(t + t) = z(t) + ri(t) + o(t 2 ). (31) 
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As we saw in Section 1.2.2, the quantity rz(<) can be much smaller than z (<). 
We can then assume that the higher order terms, in the Taylor expansion 
are negligible by comparison with the first order ones. 


As in the previous section, we denote by £ the accuracy, i.e. the solution 
is approximated with an error of the order of e. Let us assume that, on an 
interval [f, t + r], we tolerate an error of the order of £ on the approximation 
of the small scale components. From (31), we then derive an estimate on r, 
namely : 


r < 


Ke 

liWIi' 


(32) 


We denote r(e) = Ke/\ z (t) | 2 , where K is a nondimensional real constant of 
the order of the unity. As the time derivative | z (t) | 2 has a decaying spec- 
trum, the quantity t(s) then decreases with decreasing values of the level 
N{. Estimate (32) then provides a restriction on the available level on which 
the modes can be frozen, i.e. there exists a level N{ such that r(e) becomes 
smaller than the time step A t. On Figure 14, we have plotted the ratio r(e) 
for different values of jV,-. These results are obtained from the computation 
presented in Section 3. In this case, the accuracy is given by the temporal 
discretization and is of the order of Af 3 . As the time step At is of the order 
of 1 0 — 3 , we can estimate Ke ~ 10 -8 . So, for a fixed given value of the level 
A hi the corresponding time scale r(e) presents very strong variations. In a 
previous version of the algorithm, r was estimated by (vk 2 N . ) _1 and then 
was constant in time ; this choice is obviously inappropriate. Deriving an 
estimate on r(e) is essential to insure the efficiency of the algorithm. Indeed, 
we can imagine a situation where the procedure (20) and (21) allows the 
choice of a level ;V, , while the condition (32) is violated, i.e. that r is much 
smaller than the time step. Hence, the constraint (32) provides an additional 
way to determine A!,-, and N{ 2 . 


In the transition range [A r , J (<y), N i2 (tj)], we have a time estimate T Nt (e), 
given by (32), for level N t : 


T N,(e) 


Ke 


I z;v,(f) | 2 ' 

As we have noted above, the quantity r/v, (e) decreases when A, decreases. 
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According to the definition of a V-cycle, on the lower level Aq(ij), the char- 
acteristic time scale has to satisfy : 

tjv m (£)>A t. ( 33 ) 

We recall that (33) is motivated by the fact that the scales Zjv,, will be 
frozen over only one time iteration during a complete V-cycle ; (33) is then 
a constraint that N{, (tj) has to satisfy. As we have previously remarked, the 
time derivative can be evaluated by the nonlinear terms, i.e. : 

I Zyv,, b - I Qn h B( u n, u n) I2 • 

At time tj , 8 inl (y Wj) ,z Wll ) is obtained by using the following relation : 

B tn t(yN, ] ,'ZN, l ) = B(u N ,u N ) - B(y Nli ,y Nil )- 

Hence, B(un,Un) will be computed at that time. Then, the computation of 
the quantity 

/ \ & £ ^ ^ ( 34 ) 

= I (ti) b I Q^B(u N ,u N ) \ 2 

will not add extra cost. With (33), we are sure that on all levels N, higher 
than Njj , the corresponding scales can be frozen during more than one time 
iteration. 


The characteristic time r^, 2 (e) provides an estimate of the global time 
length T c (tj) of the whole cycle" [tj,tj + T c (tj)\ ; Tjv, 2 (e) can be evaluated as 
in (34), i.e. : 

T N, 2 (e) - | qN b(u n ,u n ) | 2 ‘ ^ ^ 


Remark 1 : if N i2 lies in the quasi-static range, another estimate can be 
derived ; as we have seen before, we have : 


I z/v, 2 (tj) | 2 ^ V I Az Nia (tj) | 2 • 


Then, we can write : 

I ZjV, 2 (*j) I2 > c ll" I Az I2, 


27 


where c n is a constant of the order of unity. From the definition of the norm 
| — I 2 and of z/v , we can obtain : 


" I Azyv.^j) | 2 > v ( k N , 2 ) 2 | z Nii {tj) | 2 . 


We then obtain the following estimate of T^(e) : 


TN ^ < C 11 " ( k N<2 f | z Ni3 (tj) | 2 ' 
Moreover, we recall that, by definition : 

\ZN, 2 {tj)\2= 02 (tj) | y/V, 2 (<j) | 2 • 

We finally have an a priori estimate for r/v l2 : 


t n , 2 (e) < 


Ke 

v k N, 2 2 0 2 (t } ) | y N, 2 (t 3 ) | 2 ' 


( 36 ) 


If e is the accuracy of the time scheme, we recall that K corresponds to a 
high order derivative of u /v versus time. Hence, we can reasonably assume 
that K is at least of the order of | y |: 2 • This case occurs when a Direct 
Numerical Simulation is performed. 


Time scale estimate for the transfer terms. 

Let us denote by y instead of y^r the large scale component of the flow. We 
recall that y is governed by the equation : 

y- i /Ay+Pfi(y,y) + PB int (y,z) = Pg. 

d 

Here, P denotes Pyy, and y = Let us rewrite the previous equation under 
the form : 

y = t/Ay + Pg - Pfi(y,y) - P5 int (y,z). 

We introduce the function F defined by 

F(y) = /'Ay + Pg — P P(y, y), 

so that : 

y = P(y) - PP int (y,z). (37) 
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In this form, it clearly appears that the transfer terms PB mt (y, z) is a cor- 
rection to the time derivative of the large scale components. As it was done 
before for the small structures, we derive a Taylor expansion of y(t) : 

y (f + r) = y(t) + Ty(t) + ^-y(t) + o(t 3 ). (38) 

We assume here that the terms of order larger than three are negligible. 
From (37), we can derive : 

y(0 = F( y)-PB int (y,z). 


Reporting this expression into (38), we obtain : 

2 2 

y (t + r) = y{t) + Ty(t) + jF{y)-jPB int (y,z) + o{T 3 ). 

As we tolerate an error of the order of e on y(t + t), the coupled nonlinear 
terms PB in t(y,z) can be frozen during a time r, if : 

y | PBi wt (y,z) | 2 < Ke. (39) 


Then, it follows an estimate on r : 

2Ke 


1/2 


r < 


PB int ( y,z) | 2/ 


= r'(c). 


(40) 


We note that, if £ is the accuracy of the scheme, condition (40) is necessary 
to preserve the order of the time scheme. On Figure 15, we have plotted 
the evolution of the ratio r'(e) for different levels of refinement A r ,-. As for 
r(£), the quantity r'(e) decreases when Af t - decreases, which is due to the 
fact that PBint( y, z ) has a decaying spectrum, like z(f). So, for the levels 
Ni lying in the transition range, i.e. between iV,-, and JV,- 2 , the value of r'(e) 
corresponding to the level is the most restricted one. Hence, in order to 
control the variations of FW, £?int ( y TV, , z^r t ) on the different levels, a sufficient 
condition is to estimate r\e) on the lower level Ni y of the transition range. 
We want to note that the mathematical estimates which can be derived on 
the time derivative P6 mi ( y,z) of the transfer terms do not provide efficient 
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information ; nevertheless, the numerical experiments show a correlation 
between 

| PB,,„( y,z) [ 2 | i(t) |; 

I PB- m ,(y,z) | 2 “ \z(t)\ 2 

as we can see on Figure 11. So, we deduce that 


T 

2 " 


PB int (y,z) | 2 > c 12 y | Ffi int (y,z) | 2 


z(0 b 

z(0 I 2 


where c is nondimensional constant of the order of the 
that : 


T < 


2Ks 


z(t) 


1/2 


,ci 2 | z (<) I 2 | F5 int (y,z) | 2 


unity. Hence, it follows 
= r w (e). (41) 


We then obtain with (41) an estimate of r as a function of | z (t) | 2 . From 
a computational point of view, (41) is much more efficient than (40). As it 
was noted before, we have to derive an estimate on r'L (e) in order to control 
the time variations of the transfer terms which depend on the scales in the 
transition range. We now recall that the level N it is defined by the evaluation 
of the ratio : 


01 (tj) 


II y || ' 


Moreover, this ratio is equivalent to the ratio of the coupled nonlinear terms 
of the large scales in the energy norm, namely 


I C/v,, B mt (y Nn , z Nt] ) | 2 ^ || z Ni[ (tj) || 

I flu, B(yn,,yn,) | 2 - Ca II yjv.,(C) II " 


where c is of the order of the unity. The estimates can then be written under 
the new form : 


r Ni 


(e) < 


' ‘2Ke I zjv„(C) |2 ~ 

M l i) I p N. t B( yw.^y/v,,) | 2 x | z N it (tj) | 2/ 


1/2 


(42) 


In (42), the time derivative of the small scale components Zv {tj) can be 
estimated as it was done previously. Finally, (s) provides a constraint on 
the choice of the level TV,-,, while (e) and (e) provides two estimates 
of the length of the whole cycle. 


30 



2.3 Approximate equation for the quasi static scales 
z Nl '• projection on an Approximate Inertial Man- 
ifold. 

In this section, we intend to discuss the efficiency of the Approximate Inertial 
Manifolds (AIM) and show in which part of the spectrum they can be used. 

We consider an approximation of the equation of the small scales (12) in 
which we drop the coupled nonlinear terms : 

^ - z/Az + QB(y,y) = 0. (43) 

at 

We now introduce an operator e _ " (A defined by : 

e -1/(A z = £ e"l k l 2 ‘u(M) 

ke/jv 

Then by applying this operator to (43), we obtain : 

l(e-" (A z ) = -e-^QBi y,y). (44) 

We assume, in agreement with the method described in Section 2.1, that z 
has been frozen on the time interval [tj,tj + r c (tj )] where r c (t : ) is estimated 
as in the previous section. For the sake of simplicity, we write here t instead 
of tj and t c instead of T c (tj). We then integrate (44) over the interval [t,t + T c ], 
which yields : 

Z (t + T c ) = e‘ ,TcA z(0- J t+Tc e- /{<7 - t - Tc)A QB(y(a),y(a)) d(T. 
Consider then the following approximation of the right-hand side : 
f t+T ° e _1/(<T_( ' +Tc))A QB(y(^),y(cr)) da 

Jt (45) 

~ (-i/A) _, (l - c VTcA )QB{y(t + r c ),y{t + r c )). 


Finally, z (t + t c ) is computed by : 
z(t + r c ) = e" T ‘ A z (t) 

-(-//A) _1 (l - e^ A )QB{y(t + r c ),y(t + r c )) 


(46) 
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(46) means that the small scales are slaved by the large ones. Also, from a 
computational point of view, (46) provides an efficient way to evaluate the 
small scales. In comparison to a classical time scheme, (46) presents the 
advantage that only one evaluation of the nonlinear terms is required. The 
error occuring by using (46) instead of integrating the small scales equation 
is constituted by two components, namely the time discretization on the 
approximation of the integral and the dropped terms Q B[ nX ,{y , z) : 

/ t + Tc 

e-^-^))^QB- lM (y{a),z{cT)) da | 2 

/*f+T c 

+ i J t e “*' {<,-( ‘ +Te))A (QB(y(a),y(a))~ QB(y(t + T c ),y(t + T c )))da 
— r c | QB mi (y, z) tyTc \2 + 2 r c \ QB( y,y)* jTc I2 

(47) 

where | QB mt (y y z) t}Tc I 2 = max | Q Bi nt (y(cr), z(cr)) | 2 , and similarly for 

I QB{ y,y),. n | 2 . 

Considering the previous discussions, 

I *(0 I2 ^ I QB inl (y(t),z{t)) | 2 

- | Qfiint(y,z) tjTc | 2 

- I QB(y,y)t,r c I2 • 

Then, we obtain an estimate of the error : 

I e(z) I 2 < ci 4 t c | z (t) | 2 . (48) 

Recalling that r c is defined such that 

t c | z (t) | 2 < K £, 

we then have the following estimate : 

I e(z) | 2 < c 14 I< e. 

Then, from the definition of the level jV, 2 and r c , the error introduced by 
using the approximate equation (46) is always smaller than e. We want to 
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note that the error e(z) does not depend directly on the level iV , 2 : there is 
no restriction on the value of N i2 . The Approximate Inertial Manifold (46) 
can then be used to simulate the evolution of the fine structures of the flow 
even for wave-numbers lying inside the inertial range of the spectrum. In 
the numerical simulations presented in the next subsection, (46) will be used. 

Remark 2 : Let us now consider the Approximate Inertial Manifold intro- 
duced in [1], namely : 

- t/Az + QB(y,y) = 0. (49) 

We recall that (49) consists of an approximation of the full equation of the 
small scale components z, where the time derivative z as well as the coupled 
nonlinear terms QB mt (y,z) have been dropped. In the case considered here 
of periodic boundary conditions, it is easy to invert the Stokes operator (-A) 
and then z can be evaluated as a function of y : 

z = (-i/A) _1 QB(y,y). (50) 

At this point, we note that for large values of r c , (46) and (50) are equivalent. 
The order of magnitude of the dropped terms in the small scale equations 
may induce a restriction on the use of (50) to evaluate z. In fact, we want 
to find criteria telling in which range of the spectrum (50) can be applied. 
The spatial error S( z) appearing when (50) is used is exactly given by the 
difference between (50) and the z equation, i.e. : 

S( z) = (i/A) _1 (z + QB int (y,z)). 

We can then estimate 

I 6{z) | 2 < {vk 2 Ni ) _1 | z + QB int (y, z) | 3 . 

In the numerical experiments that we have conducted, we have seen that the 
right-hand side of the previous inequality is of the order of (t^^,) 1 | z | 2 ; 
hence it follows that 


| s(z) | 2 < C15 (j'A:)^) 1 | z(f) | 2 , 

where Cjs is a nondimensional constant of the order of unity. Then, the spatial 
error introduced by (50) mainly depends on the size of | z (f) | 2 . Assuming 
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that Ni lies in the dissipation range and that | i(t) | 2 can be estimated by 
| /'Az | 2 , we obtain : 

I <H Z ) I2 < c 16 | z | 2 , 


by using 

| i Nl 1 2 ~ (i'k 2 Ni ) | z Nl | 2 . 

Hence, in that specific case, the error S(z) is of the order of z itself. Thus (50) 
can be used in the quasi-static range of the spectrum where z is of the 
order of the scheme accuracy. On Figure 16, we can see that the quantity 
I z(0 b becomes larger than | z | 2 itself when the cut-off value Ah 
decreases. Hence, it seems that if (50) is an efficient way to compute the 
very fine structures of the flow lying far inside the dissipation range, it is no 
longer the case for the scales of the order of, and immediately larger than 
the dissipation scale £ v . 


2.4 Description of the complete multilevel algorithm. 

In this section, we summarize the complete multilevel method which includes 
the time scales derived in Section 1.2.2. We still denote by e the accuracy 
of the computation ; we recall that e is a given parameter in the following 
algorithm. 

As in subsection 2.1, we choose a sequence of levels Ni such that : 

< N 2 < ... < Ni < N i+ ] < ... < N. 

The whole time interval of the simulation, namely [f 0 ,<o + T], is split into 
several intervals of the form 


[tj, tj + , 


j-1 

where tj = t 0 + T c (t k ). Futhermore, we assume that the final time t 0 + T 
k= 1 

satisfies 

t 0 + T = £ m , 

so that we have in intervals. Let us now assume that the approximate solution 
U;v(x,£) is known at time tj with j < m. As in 2.1, we compute two levels 
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of refinement N.^ij) and N l2 (tj) according to the procedures (20) and (21). 
Moreover, we impose that 


T N tl ( £ ) 


K £ 

ZJV M b 


> A*. 


With (35) and (42), we derive an a priori estimate of respectively T Ni 2 (e) and 
t'x (e) ; we then obtain an evaluation of the length r c (/,j) of the j— cycle : 

T c (tj) = min(T N<2 (£), r^ tj (e)). (51) 

We note here that, with this definition, T c (tj ) can be smaller than one V- 
cycle, i.e. ( 2 (i 2 -ii) + l)A* ; in such case, levels N it and N h are too small and 
need to be reevaluated. With (51), we have an a priori estimate of the global 
length of the whole cycle. Finally, we have computed the three characteristic 
values : 

N i2 (tj) and T c (tj). 

As in 2.1, the integration is performed on the interval [tj,tj + T c (tj)] by a 
succession of V— cycles. At the end of each V— cycle, i.e. at time tj+pV 
t- 4 - (2 p(i 2 — *i) + 1 )A£, we derive an a posteriori estimate of the quantities 
r" (£) and r^(e). Hence, we take into account the evolution of the scales 
lying in the transition range of the spectrum of the velocity (see Figure 3). At 
this time, if {ij+T c {tj))-tj+ p V is larger than one full V-cycle, i.e. 2 (* 2 — )~F 1 

time iterations, we perform another V-cycle after reajdusting the value of 
Tc (tj) with the new values of (e) and r Nia (e). Now, in the other case, i.e. 
if (tj + T c (tj)) — tj+pv is smaller than one V-cycle, we stop the whole cycle 
by saying that 

tj T c (tj ) tj-\- 1 - 

At tj + T c (tj), we compute the small scale components z/vi 2 of the spectrum 
by projecting the solution on the Approximate Inertial Manifold (46). 

Then, we readjust the two levels and restart a new cycle as it was done in 
Section 2.1. 


Before we conclude this section, we want to note that with this algorithm, 
in opposition with the previous ones (see for instance [10], [7] or [17]), the 
constants 9\ and 0 2 of the procedures (20) and (21) can be fixed very easily. 
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Indeed, as we have previously seen the parameter 9 2 is chosen so that e( Zjv, ) 
is of the order of e 2 . The constant which provides an estimate of the ratio 
II z/v,, || / || y n h II, can be evaluated at the initial state, i.e. t = 0, by 


At lly^ III Yn h | 2 ' 

Hence, we insure that : 

A/v m = At | z Nii | 2 < cAt | yjv M | 2 || z Nij ||< e. 

Then, the condition on r N) ^ e ) is satisfied i.e. t Nii ( e ) > At. Moreover, we 
can implement a self-adaptative procedure allowing a dynamical reevaluation 
of these constants during the time evolution. So, if 9 , and 9 2 were previously 
fixed in an empirical way in the algorithms, we have found now a more 
efficient way to evaluate these constants. 
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Figure 4: Time evolution of the ratio | ZjV l r for Ah — 32, 64, 128 and 
6 I YNi |2 

Ah = 196. 


0.159E+01 


0.536E-02 J\ 




0.181E-04 

» \ /*'“ ii- /* >, 

! \ /' W '-..'jj-S-j *"V _ 

\ Vw-^ \ /wyw u A‘ 

0.610E-07 . Vl . r -x 

i r v V A \ 


\ J 

0.206E-09 . \ jT 


V ^ 

\/Vr- J 

i/ 


0.694E-12 


20 40 60 80 100 


Log. Scale 


Figure 5: Time evolution of the ratio '| ^ Nl 

32, 64, 128 and N l = 196. 
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Figure 6: Time evolution of 
yVi = 196. 
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Figure 7: Time evolution of | z Nl | 2 and | Q% x B int ( y Nl ,z Nl ) | 2 , for /V, = 32 
and Aj = 64. 
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Figure 8: Time evolution of A t \ z^ l I 2 and | z^v, I 2 for — 32,64, 128 and 


196. 
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Figure 12: Time evolution of | z Ni | 2 , // | Az Nl | 2 , | Q^B{ y N] ,y Wl ) | 2 and 
I | 2 for /V, = 32, 64, 128 and 196. 
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Figure 16: Time evolution of (t'fcjv,) 1 I Zyv, I 2 and I zjv, I 2 for M = 
32, 64, 128, and 196. 
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Numerical results. 


3.1 Comparison with the Galerkin method and with 
a previous version of the multilevel method. 

In this section, we report numerical results obtained by using the nonlinear 
Galerkin method and a classical Galerkin projection. The flow of Kolmogorov 
type is forced by a time independent external force f which acts, in the 
spectral space, on only some low frequency components of the velocity field. 
The initial condition is chosen so that its spectrum has a specified shape but 
the phases of its Fourier components are randomly chosen. So, the flow at 
time t = 0 has no organized struture. We have let the flow evolve on over 
10 5 time iterations, i.e. from t = 0 to t = 100 ; this is much longer than the 
integral time scale, which is of the order of the unity. We have compared the 
solutions obtained with both methods. 


3.1.1 Description of the computation. 

The initial condition here is computed from a given spectrum of the initial 
vorticity — V x u^, where is the following expansion : 


u o V ( x ) = £ u 0 k (i)e* kx , 

k 6 /jv 

with x = (aq, £ 2 ) € ST We choose u>q by setting 

^o,k = | «>o,k I e ,e k, 


(52) 


(53) 


where 0 k £ [0,27 t] is generated by a random function, and : 


^o,k 


(k + ( V ^) 5 ) 1 ' 2 
k 0 


if k =| k |< k a , 


otherwise ; 


(54) 


c \ 7 is determined such that | Lotf |l°°= 2.0 ; k a is equal to 60. At this point, 
we note that if | u) o k |~ ck^ , then the energy spectrum of Uo 

Eo(k) — y; (I Uo,k | 2 + | ^ 0 ,k | 2 ) 

Ikl = k 
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is like ck 2(3 ~ 2 . So, in the case presented here ft = —0.5 and we have Eo{k) ~ 
k~ 3 . On Figure 28, one can see the isovorticity lines of the initial velocity 
field and Figure 17 shows the energy spectrum E(k). 

The external force f is constant in time and has only a few non-zero 
wavenumbers, namely : 

fk = (/i,k?/2,k)i 

with 

r i ; 1>k i = c 18 if k e z v \k,\ + \k 2 \= 3, 

( | J lk | = 0 otherwise, 

d 8 is determined such that | f | 2 = 0.225. The Fourier coefficients of f are 
finally obtained by 

fi , k = I Ak I e"k, 

where the phases 0 k € [0, 27 t] are randomly generated. 

In order to describe all the scales of motion, the number of modes N 
in each dimension of the space must be chosen so that the associated grid 
size 2 n/N is smaller than the dissipative (Kraichnan) scales ( v ; in term 
of wavenumber, it means that k n > k v . We recall that under the dissipative 
scales, the motion is damped by viscosity. In fact, the total number of degrees 
of freedom needed to describe the motion, from the dissipative scales to the 
large scales containing eddies, can be estimated by the ratio ( k v /k 0 ) (see [22] 
and [16]). Constantin, Foias, Manley and Temam, in [16], have related this 
quantity to the dimension of the attrator of the Navier-Stokes equations 

(k v /k L ) 2 ~ Rcl, 

where Rei is the integral scale Reynolds number, which can be defined by 



v 


Here t? 2 = (2/ | fi |) e(u) and L = l/k^is the integral length scale defined 

as : 

Y 1/2 

L = — 7TT, where E = (1/2) d 2 and t) = // | Au | 2 / | fi | (t^ = 10 ). 

7/ 1/J 
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L is of the order of 1.72 and t? is of the order of 0.45 (see Figure 19). It 
follows that Re i = 784 and k v = 17, which corresponds to N v = 2 k v = 34 
modes in each direction of space. In order to be sure to resolve all the scales 
of motion, we have chosen N = 256. Figures 17 show the energy spectrum 

m : 

E(k,t) = £ |u(k,()| 2 , 

|k| = k 

for various times. It seems that the dissipation wavenumber k t) is of the 
order of 20. Hence, the previous estimate based on the dimension of the at- 
tractor matches well the computational results. We note also that, if the 
phenomenological theory of turbulence of Kraichnan, predicts a decay of the 
energy spectrum like k~ 3 , results presented here seem to show a faster decay 
closer to k~ 4 . This results are in agreement with the ones obtained by Orszag 
in [23] and by Brachet et al. in [24], which show that a k~ 3 energy spectrum 
can only be obtained when the Reynolds number is larger than 25, 000. 

We also remark, that if, at t = 0 the small scales corresponding to a wavenum- 
ber larger than 60 are set to zero, a dissipation range appears very quickly, 
as we can see at t = 5. The enstrophy transfer from the large scales to the 
small ones acts on the smallest scales after a few iterations. Figure 18 indi- 
cates that the transition period is very short. The small scales are damped 
by viscous effect until an equilibrium between viscous and nonlinear terms 
appears. After that, we can see that | Zjv, I 2 oscillates in time and seems to 
become completely independent of its initial value. Figures 28 and 30 repre- 
sent the isolines of the vorticity, at different times in the interval [0, 100]. We 
can see that the very small random structures of the flow at the initial time 
disappear quickly. It appears that fusions of these very small structures lead 
to larger ones. So, after a transient period, the flow is mainly constituted by 
large structures. 

The time step is chosen by considering the accuracy and the stability of 
the computation. For the stability, At must satisfy a CFL condition like 

At N I Uyv | L oo < a (< 1). (55) 

From Figure 25, we can see that | ujv |l°°< 1-25 during all the computation, 
so (55) implies the following restriction on the time step : 

At < 1.7 1(T 3 (a = 0.5). 


48 



We have set At to 10 -3 . Here, the smallest scales (Figure 18) are of the order 
of lO -10 . Hence a time step of 10 -3 allows to recover most of the spectrum. 
Indeed, the time differentiation scheme used here is a third order method, 
then the accuracy is of the order of 10 -9 . 

3.1.2 Comparison with a previous version of the algorithm. 

In a first time, we present results obtained with a previous version of the 
multiscale method. In this version of the algorithm, r c was set to (i/ ^iv 12 ) 
and the constants d u 0 2 of procedures (20) and (21) were a priori chosen 
(see [8]). This multiscale method is compared with the pseudo-spectral 
Galerkin method. 

To perform this analysis, we have retained two different points of the compu- 
tational domain, namely x\ — x 2 = jy- (-y — l) a °d x-i — x 2 = ^ (-5 l), 

and we have stored the value of the horizontal component of the velocity 
ui(xi,x 2 ,t) at these points during the time evolution. On Figures 20 and 21, 
we have plotted the time history over the interval [0, 100] of those two char- 
acteristic values of the flow. Results plotted here seem to be identical for 
both methods, but the differences between the trajectories obtained with the 
different algorithms are too small to appear on such graphic representation. 
So, we have listed below the exact values at different intermediate times for 
these trajectories. 

In Tables 1 and 2, we have listed the values of the horizontal component 
of the velocity at times t = 0,25,50,75 and 100. These results correspond 
to the Galerkin method for the first column. In the third column, we can 
see the difference between both orbits. It appears that this quantity grows 
as time evolves and becomes much larger than the accuracy which is of the 
order of At 3 = 10~ 9 for this computation. On this computation, the current 
level Ni(t) oscillates between 108 and 200. Looking backward to the results 
presented in section 2.2 on the estimates of the characteristic length t c (see 
Figures 22), we note that : 

| t N)j < 10 -3 = At, 

( < 10 -2 = lOAt, for N<, = 108. 

Hence, levels JV ,-j used by the multiscale method are not appropriate. 
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Table 1 

Galerkin 

Multiscale 

Difference 



Version 1 


t = 0 

-0.9650082490 10" 2 


0.000 

t = 25 

0.2980770146 

0.2980770543 

3.88 10~ 8 

t = 50 

0.3592114126 

0.3592114913 

7.87 10- 8 

t = 75 

0.3819454889 

0.3819452151 

2.738 10" 7 

t = 100 

—0.123374601 1 

-0.1233715822 

3.0189 10~ 6 


Table 2 

Galerkin 

Multiscale 

Difference 



Version 1 


t = 0 

0.4501 442245 10- 1 


0.000 

t = 25 

0.3879730195 

0.3879729977 

2.18 10“ 8 

t = 50 

0.2458823192 

0.2458824078 

9.114 10~ 7 

t = 75 

-0.3332440190 

-0.3332451088 

1.091 10" 6 

t = 100 

0.3196664945 

0.3196686370 

2.142 10~ 6 


Indeed, we recall that Tty, < At means that, for the level N tl equal to 
108, the scales smaller than t tv, , can not be fixed even on one time iteration. 
AlsotyV., < 10 At means that the coupled nonlinear terms P N>1 B(y NtJ , z/y, ,) 
can not be frozen on a time interval longer than lOAt. Recalling now that 
in this version of the multiscale algorithm, the characteristic length was es- 
timated by 1 1 therefore the length of the frozen period oscillated 

between 50 At and 80 At. The constraints on the levels TV,, and /V,- 2 imposed 
by the different estimates derived in section (2.2) are violated in this compu- 
tation and so, the expected accuracy can not be recovered by the multiscale 
method. 


Table 3 

Galerkin 

Multiscale 

Difference 



Version 2 


t = 0 

-0.9650082490 10~ 2 


0.0000 

t = 25 

0.2980770146 

0.2980770146 

< 1.0 io- 10 

t = 50 

0.3592114126 

0.3592114126 

r< l.o io~ 10 

t = 75 

0.3819454889 

0.3819454893 

4.0 10~ 10 

t = 100 

-0.1233746011 

-0.1233746004 

o 

h— * 

o 

1 

o 


50 






Table 4 

Galerkin 

Multiscale 

Difference 



Version 2 


t = 0 

0.4501 442245 10 _1 


0.000 

t = 25 

0.3879730195 

0.3879729977 

2.0 10- 1U 

t = 50 

0.2458823192 

0.2458824078 

< 1.0 io~ 10 

t = 75 

-0.3332440190 

-0.3332451088 

2.0 10" 10 

t = 100 

0.3196664945 

0.3196686370 

< 1.0 10" 10 


3.1.3 Comparison with the improved version of the algorithm. 

In Tables 3 and 4, we report results similar to those in Tables 1 and 2, 
but now, the second column corresponds to the version of the multiscale 
algorithm presented in Section 2.4. Here, the trajectories of the multiscale 
method remain close to the trajectories obtained with the classical method. 
The difference is less than A t 3 = 1(T 9 over the whole time interval [0, 100]. 
Here, the level Nn is always larger than 128, so that the estimated limit 
value of t c is greater than At (see Figures 22 and 25). The levels JV tl and N l2 
chosen by the new algorithm are higher than those obtained by the previous 
version. This fact is due to the restrictions imposed by the new criteria on 
t c . For some values of the time t, the level Afj, decreases to a lower level 
for a short time interval and then goes up to its last value. In such a case, 
we have remarked that the restriction imposed on r c induces a restriction on 
the level N iv Indeed, even if a level is acceptable on a few multigrid cycles, 
variations of the smallest scales or of the transfers terms become too large, 
so that N { i has to change to an upper level. Moreover, we want to mention 
that, during the whole computation, neither the restriction due to the small 
scales evolution t Ni , nor the restriction due to the transfer terms evolution 
t" m is dominant. Therefore, it is necessary that both criteria be retained. 
Finally, we want to note that the algorithm can still be improved. Indeed, we 
recall that the actual choice of levels TV,, and N h are determined according 
two different criteria. One is based on the estimate of ratios of the kinetic 
energy (or enstrophy) of the small scales over the kinetic energy (or enstro- 
phy) of the large ones, and the other one is based on the estimates of the 
critical characteristic time r c . So, it follows that r c may have a relatively 
large value, and then, the levels should be adjusted to a lower value in order 
to have a more reasonable estimate of t c . This can be viewed on Figure 27, 
where the time evolution of r c is plotted. In fact, the very strong oscillations 
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of t c are not always necessary, and, by optimising the evaluation of /V u and 
/V, 2 , this time evolution can become smoother. 

Finally, to achieve the comparison between the different algorithms, we want 
to note the non-negligible fact that the multiscale method requires twice 
less CPU time than the classical pseudospectral method. On Figure 24, the 
quantity : 

T NL g - To 
To 

is plotted, where Tnlg is the CPU time required by the nonlinear Galerkin 
method (multiscale), and Tq the CPU time required by the classical Galerkin 
method. 

Finally, we want to mention that these simulations required more than 75 
hours of CPU on a Cray2, without counting all the preliminary tests needed 
to the developments and improvements of the algorithm. 
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Figure 17: Energy spectrum at t — 0, 5, 30, 50, 85, 100. 
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Figure 18: Time evolution of | Zjv, | 2 for JVj = 32, 64, 128 and TV, = 196. 
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Figure 19: Time evolution of j u^(t) | 2 . 
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Figure 20: Time evolution of the horizontal component of the velocity 
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Figure 21: Time evolution of the horizontal component of the velocity 

u N (x,y,t) at point x = y = - !)• 
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Figure 24: Time evolution of the CPU time for the classical method (full 
line) and for the multilevel method (dashed line). 
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Figure 25: Time evolution of the levels | u;v(t) | L „ 
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Figure 27: Time evolution of the characteristic time T c (tj) 
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Figure 28: Vorticity structures at time t — 0. 
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Figure 29: Vorticity structures at time t = 30. 
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Figure 30: Vorticity structures at time t = 60. 
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Figure 31: Vorticity structures at time t = 100. 


1.7779 

0.8836 

- 0.0107 

- 0.9050 

- 1.7994 



60 



3.2 A direct numerical simulation at higher Reynolds 
number 

We present here numerical results obtained in a numerical simulation, sim- 
ilar to the previous one, but with a larger value of the Reynolds number. 
The external force is kept the same and the viscosity is divided by 4. The 
initial condition is still a random field computed from a given vorticity. The 
computation is performed over 50,000 iterations, so that, the flow is no more 
dependant of its unstructured initial state. An a posteriori analysis of the 
behavior of the adaptative multilevel procedure confirms the previous as- 
sumptions and proves that the multiscale method is well adapted to Direct 
Numerical Simulations. This computation required about 50 CPU hours on 
a CRAY2. This total computing time includes all post-treatments done by 
the code, i.e. computations of different norms of several quantities related to 
the small and large scales. The CPU time spent to compute the velocity field 
is 32 hours, which corresponds to 7 10 6 second per mode and per iteration. 


3.2.1 Description of the initial condition. 

As in (3.1), the initial field u 0 is computed in the spectral space from the 
coefficients of a given vorticity 

ug = £ u> 0 (k)e ik ' x , (56) 

keZ 2 ,|k|</v 


and the coefficients w 0 (k) are given by 


u;o(k) 


c 19 | k |-5 e i0k if | k | < k a 
0 otherwise 
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(57) 


where c 19 is such that | \ LOO = 2.0. As we can see on Figure 32, the slope 
of the energy spectrum at the initial time is equal to —3. 

The viscosity is set to 2.5 1(T 4 , which gives a Reynolds number equal to 

6,328 and the integral scale L = 2.93. So, the dissipative wavenumber k v is 
of the order of 28 which corresponds to N n = 56. As we want to perform a 
Direct Numerical Simulation a total number of modes of 512 in each direction 
provides a grid fine enough to resolve the scales under the dissipative ones. 
On Figure 32, we can see that the dissipative wavenumber obtained by the 
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computation is below 60. Moreover, on Figure 32 we see that the slope of the 
energy spectrum is smaller than -3, which is a faster decay of the energy 
than for a fully develloped turbulent flow. 

As in the previous computation, a time step equal to 10 -3 is small enough to 
insure the stability of the scheme and allows to compute all the scales with 
enough accuracy. 

3.2.2 Analysis of the computation. 

We first want to make some remarks on the behavior of the multiscale 
method. Figure 36 shows the time evolution of the two characteristic lev- 
els Ni x and N, 2 , which define the transition range. As it was expected, the 
variations of N i} and N h follow the variations of the ratios 

I b and I flv,£My;v,,Zjv,) | 2 
I 1 2 I PN,B(y Nl ,y Nl ) | 2 

for values of N\ larger than 256, as we can see on Figures 33 and 34. On 
Figures 37 and 38, the variations of the characteristic times t/v, and r ^ are 
plotted. We remind that Tty is used to determine the lower level and the 
length of the period during which the smallest scales, i.e. for t < l N . , can be 
frozen. From figure 36, we note that the lower level of the transition range, 
namely N it , has to be larger than 256. In fact, we find N i{ of the order of 
320, for its lowest value. Figure 35 confirms that the time derivative of the 
velocity has a decaying spectrum. 

On Figures 39, we have plotted the time evolution of the different terms 
appearing in the equation of the small scale components Zjv,. As we have 
observed in the previous computations, the time derivative | zjv, I 2 is of 
the order of the coupled nonlinear terms | Q * B tnt (y Nl , ZyV| ) | 2 , while the 
dissipative norm v | Az/v, | 2 is much smaller for lower values of N\. For 
the scales lying far inside the dissipation range, we also note that all the 
quantities are of the same order. 

Figure 40 show the vorticity at different time of the integration interval [0, 50]. 
The unstructured initial field disappear completely after a fiew times the eddy 
turnover time. After a short transient period, the flow evolves by keeping 
the same global structures. Figure 40 represents three-dimensional views of 
two-dimensional maps. Lighting technics have been used, so that shadow 
effects allow to see the very fine structures of the flow. 
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Figure 32: Energy spectrum at t = 0, 15, 30, 50. 
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Figure 39: Time evolution of (resp.) v | Azjv, I 2 , | zyv, h, I Q^ x B{y^^yN x ) I 2 
, and | Q^BU y Nl ,z Nl ) | 2 for N, = 64, 128, 256 and N x = 480. 



Figure 40: Vorticity structures at time t = 30. 


fc 


- 3.6502 
| 1 . 7402 


i 


- 0.1697 

- 2.0797 

- 3.9897 



Figure 41: Vorticity structures at time t = 50. 
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3.3 Simulation of the whole dissipation range with 
the multiscale method. 

In this section, the example is the same as in the Section 3.1. Our purpose 
is to study the effect of the cut-off value TV,, on the computed solution, so 
that several numerical experiments have been conducted. We have decreased 
JV M so that it is of the order of the dissipation wavenumber k v ; the whole 
dissipation range is then modelized by the multiscale strategy described in 
previous sections. In both simulations, we have noted that the multilevel 
method allows to recover all the large structures of the flow. Finally, we have 
made a similar test with the classical method, i.e. we have decreased N in 
order to estimate the lower level required to recover the large scales. 


3.3.1 Analysis of the numerical simulation. 

In order to decrease the value of N tl , the parameter e has been set to 10 
In this case, the level ZV„ adjusts itself, using the procedure described in Sec- 
tion 2, to the value 24 which is smaller than N v = 34. Figure 43 shows the 
evolution of the two characteristic levels and fV,- 2 , and Figure 44 shows 
the evolution of the characteristic time r c . The level N n is quasiconstant and 
is equal to 24 while the level N t2 is approximately equal to 128 ; iV, 2 is chosen 
to be far inside the dissipation range. On Figures 42, we have represented 
the energy spectrum of the computed solution at different intermediate times 
of the interval [0, 65] on which is conducted the computation. There is no 
energy pile-up at high wavenumbers and the dissipation of the enstrophy is 
well preserved. By comparing the vorticities obtained with this simulation 
and the Direct Numerical Simulation performed with the classical Galerkin 
(pseudospectral) method with a spatial resolution equal to (256) 2 , we note 
that all the large scale vortices are well described by the multilevel method. 
At the times t = 60 and t = 65, we still have the fusion previously mentioned. 
We also note that oscillations appear on the domain. They are due to the 
approximation made on the small scales. They also point out the problem 
of separation of scales related with a Fourier decomposition. Indeed, Fourier 
waves oscillate over the whole domain and they can not be directly associ- 
ated with a scale vorticity. This nonlocal property implies the oscillation 
appearing on Figures 47. Nevertheless, it appears clearly that all the large 
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Table 7 

u[ !, ’(x u x 2) 

u i NL (x\,- r 2 ) 

Difference 

t = 10 

0.8299845080 10“’ 

0.8355520907 lO -1 

5.56 10- 4 

t = 20 

0.2091267106 

0.2087057712 

4.21 10- 4 

t = 30 

0.3664263833 

0.3652791132 

1.14 10- 3 

t = 40 

0.2926552919 TfF 7 

0.2587214575 10“' 

3.39 10- 3 

II 

Oi 

o 

0.2458823192 

0.2432537325 

2.63 10" 3 

t = 60 

0.1573811778 

0.1414418124 

1.59 10~ 2 


Table 8 

uF''(x,,x. 2 ) 

u( \^ L {x^■,X2) 

Difference 

t = 10 

-0.7118241214 10- 1 

-0.7088992727 10" 1 

2.92 10~ 4 

t = 20 

0.1478691156 

0.1470191181 

8.49 10“ 4 

t = 30 

0.4236670791 

0.4264640060 1 

2.79 10- 3 

t = 40 

0.4898684126 

0.4800345247 

9.83 10~ 3 

t = 50 

0.35921 14126 

0.3480702322 

1.11 10-' 2 

/ = 60 

0.5173406258 

0.5179894219 

6.48 10~ 4 


structures are captured with the multilevel method. 

The total CPU time required for this simulation is about 4796 seconds and 
the difference in l? (or energy) norm with the solution obtained with the 
previous DNS is of the order of .16 10~ 3 at t = 50. The upper level l\ of 
the multilevel procedure is approximatively equal to 128 during the whole 
computation. The small scale energy on this level is less than 10 -8 and then 
much less than c. Hence, the level N l2 can be decreased to at least 64 as we 
can see on Figure 16 ; in that case, the CPU time will be reduced. 

On tables 7 and 8, we have compared the solutions obtained here with 
the one obtained by direct simulation. It appears that the difference between 
the two solutions is of the order of 10~ 2 and hence, of the order of f. The 
oscillations, appearing on figures 45, 46, and 47 are probably due to the fact 
that the intermediate scales, lying in the inertial range and corresponding 
to wavenumbers larger than , are not relaxed at each time iteration. A 
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better understanding of the effects of the multilevel procedure on these scales 
must be done in order to improve the results presented here. 

3.3.2 Comparison with a low resolution Galerkin simulation. 


Finally, in order to show the efficiency of the multilevel method on such sim- 
ulations, we have performed an additional test with the classical method. 
Indeed, we have first fixed the total number of modes to 32 and we have 
integrated the system over 50,000 time iterations, which corresponds to the 
time interval [0,50]. The solution obtained is different than the one obtainec 
when a Direct Simulation is done, i.e. with 256 modes i Figure 3 
the vorticity structures obtained with this simulation at time t - oO. llie 
number of modes is not sufficiently large so that the large scales can not 
be computed. The vorticity looks like a non-structured quan i T- ,g ’^ 
shows the energy spectrum at different time of the interval [0,o0], By fix- 
ing the resolution to (32) 2 , it seems that we do not allow the appeal an 
small scales and their dissipation mechanism occur.ng in the viscous range, 
the dynamic of the flow is drastically modified. We note that this simu ation 
is stable in the sense that no quantity grows artificially. A similar simulation 
lias been done with N = 48 instead of N = 32. The picture of the flow (se 
Figures 47, 49) is more realistic than in the previous case. The large scales 
of the flow are almost well captured. Nevertheless, the simulation performed 
with the multilevel method with a lower level N t , = 24 gives a better quah- 
tative result. The CPU time required by this simulation is 2730 sexonds ai 
the difference with the DNS is of the order of 3.1 10“* in L 2 norm. Hence, the 
accuracy recovered here is the same as the one obtained with the multilevel 
level method when N u = 24. The result obtained here confirms the fact that 
in the previous computation, the level N X2 can be decreased to 48 ; in such 
a case, the CPU time required by the multilevel method will be at least two 

times less than 2730 seconds. 


To conclude, we have proved with these experiences that the multilevel 
method allows to model ize the whole dissipation range an ie en 
inertial range without disturbing the large scales of motion. Moreover, we 
have shown that if the small scales lying in the dissipation range can e 
modelized, as well as their interaction with the large scales, they are essential 
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to describe the evolution of the flow. 
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1 = 20.0 


Figure 42: Energy spectrum at t — 5, 20, 35, and 50. 
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Figure 45: Comparison 


of the vorticity structures at time t — 60. 
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Figure 46: Comparison of the vorticity structures at time t = 65. 
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Figure 47: Comparison of the vorticity structures obtained with the classical 
flalerkin method for N = 256 and N = 32, at time t o . 
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Figure 48: Energy spectrum at t - 5, 20, 35, and 50. 
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Figure 49: Comparison of the vorticity structures, between results obtained 
with the classical method for N = 256 and N = 48, at time t = 50. 
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Figure 50: Energy spectrum at t — 5, 20, 35, and 50 
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Conclusions 


In this article, we have proposed a multilevel method which treats differently 
the large and the small scales of homogeneous isotropic flows. Moreover, 
we have derived mathematical estimates for all the parameters (cut-off level, 
time scales) involved in this dynamical procedure leading to a completely 
self-adaptive procedure. Firstly, several computations have been conducted 
in the context of DNS, i.e. the whole spectrum up to the dissipative scale was 
simulated. In such case, the multi-level method is able to recover a velocity 
field with the same spatial resolution than the Galerkin method, but with a 
substantial speed-up of at least 2 in CPU time. 

Secondly, in an approach similar to LES, we have decreased the number 
of modes retained for the resolved scales and used the same algorithm to 
modelize the interaction between the low and high frequencies. In such case, 
we have seen that when the resolved scales are larger than the dissipative 
ones, but of the same order of magnitude (1/12 as compared to 1/17), the 
large eddies of the flow are captured correctly; this is not the case when the 
pseudo-spectral Galerkin method is used with only 32 2 modes, i.e. with a 
larger scale of the order of 1/16. 

Although, we only presented computational for moderately large Reynolds 
numbers, where the small scales are strongly dependent on the energy- 
containing eddies, we think that these results are very promising. Indeed, 
we can reasonably hope that the multilevel (nonlinear Galerkin) method can 
be used efficiently for Large Eddy Simulations. Results on this point will be 
presented elsewhere, in the case of two and three dimensional homogeneous 
isotropic flows. 

Finally, we want to observe that, by opposition with LES methods, no as- 
sumptions on the energy spectra or on the velocity correlations are made here. 
Therefore the method can be applied to the simulation of non-homogeneous 
flows. First attempts have been made in this direction, for the channel flow 
problem. Our efforts will be concentrated on this problem in the near future. 
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